home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Software Vault: The Gold Collection
/
Software Vault - The Gold Collection (American Databankers) (1993).ISO
/
cdr53
/
acmalg01.zip
/
ACM526.FOR
< prev
next >
Wrap
Text File
|
1993-01-01
|
61KB
|
1,799 lines
C ALGORITHM 526, COLLECTED ALGORITHMS FROM ACM. THIS WORK
C PUBLISHED IN TRANS. MATH. SOFTWARE, 4(2), PP. 160-164.
C PROGRAM TTIDBS(OUTPUT,TAPE6=OUTPUT)
C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO-
C GRAM PACKAGE. ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS
C ARE EXPECTED TO BE ZERO.
C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS
C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C DECLARATION STATEMENTS
DIMENSION XD(30),YD(30),ZD(30),
1 XI(6),YI(5),ZI(6,5),
2 ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5),
3 IWK(1030),WK(240)
DATA NCP/4/
DATA NDP/30/
DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6),
1 XD(7), XD(8), XD(9), XD(10),XD(11),XD(12),
2 XD(13),XD(14),XD(15),XD(16),XD(17),XD(18),
3 XD(19),XD(20),XD(21),XD(22),XD(23),XD(24),
4 XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/
5 11.16, 24.20, 19.85, 10.35, 19.72, 0.00,
6 20.87, 19.99, 10.28, 4.51, 0.00, 16.70,
7 6.08, 25.00, 14.90, 0.00, 9.66, 5.22,
8 11.77, 15.10, 25.00, 25.00, 14.59, 15.20,
9 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/
DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6),
1 YD(7), YD(8), YD(9), YD(10),YD(11),YD(12),
2 YD(13),YD(14),YD(15),YD(16),YD(17),YD(18),
3 YD(19),YD(20),YD(21),YD(22),YD(23),YD(24),
4 YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/
5 1.24, 16.23, 10.72, 4.11, 1.39, 20.00,
6 20.00, 4.62, 15.16, 20.00, 4.48, 19.65,
7 4.58, 11.87, 3.12, 0.00, 20.00, 14.66,
8 10.47, 17.19, 3.87, 0.00, 8.71, 0.00,
9 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/
DATA ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6),
1 ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12),
2 ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18),
3 ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24),
4 ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/
5 22.15, 2.83, 7.97, 22.33, 16.83, 34.60,
6 5.74, 14.72, 21.59, 15.61, 61.77, 6.31,
7 35.74, 4.40, 21.70, 58.20, 4.73, 40.36,
8 13.62, 12.57, 8.74, 12.00, 14.81, 21.60,
9 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/
DATA NXI/6/, NYI/5/
DATA XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/
1 0.00, 5.00, 10.00, 15.00, 20.00, 25.00/
DATA YI(1), YI(2), YI(3), YI(4), YI(5)/
1 0.00, 5.00, 10.00, 15.00, 20.00/
DATA ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1),
1 ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2),
2 ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3),
3 ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4),
4 ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/
5 58.20, 39.55, 26.90, 21.71, 17.68, 12.00,
6 61.58, 39.39, 22.04, 21.29, 14.36, 8.04,
7 59.18, 27.39, 16.78, 13.25, 8.59, 5.36,
8 52.82, 40.27, 22.76, 16.61, 7.40, 2.88,
9 34.60, 14.05, 4.12, 3.17, 6.31, 0.60/
DATA LUN/6/
C CALCULATION
10 MD=1
DO 12 IYI=1,NYI
DO 11 IXI=1,NXI
IF(IXI.NE.1.OR.IYI.NE.1) MD=2
CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI),
1 ZI1(IXI,IYI),IWK,WK)
11 CONTINUE
12 CONTINUE
15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK)
DO 17 IYI=1,NYI
DO 16 IXI=1,NXI
DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI))
DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI))
16 CONTINUE
17 CONTINUE
C PRINTING OF INPUT DATA
20 WRITE (LUN,2020) NDP
DO 23 IDP=1,NDP
IF(MOD(IDP,5).EQ.1) WRITE (LUN,2021)
WRITE (LUN,2022) IDP,XD(IDP),YD(IDP),ZD(IDP)
23 CONTINUE
C PRINTING OF OUTPUT RESULTS
30 WRITE (LUN,2030)
WRITE (LUN,2031) YI
DO 33 IXI=1,NXI
WRITE (LUN,2032) XI(IXI),(ZI1(IXI,IYI),IYI=1,NYI)
33 CONTINUE
40 WRITE (LUN,2040)
WRITE (LUN,2031) YI
DO 43 IXI=1,NXI
WRITE (LUN,2032) XI(IXI),(DZI1(IXI,IYI),IYI=1,NYI)
43 CONTINUE
50 WRITE (LUN,2050)
WRITE (LUN,2031) YI
DO 53 IXI=1,NXI
WRITE (LUN,2032) XI(IXI),(ZI2(IXI,IYI),IYI=1,NYI)
53 CONTINUE
60 WRITE (LUN,2060)
WRITE (LUN,2031) YI
DO 63 IXI=1,NXI
WRITE (LUN,2032) XI(IXI),(DZI2(IXI,IYI),IYI=1,NYI)
63 CONTINUE
STOP
C FORMAT STATEMENTS
2020 FORMAT(1H1,6HTTIDBS/////3X,10HINPUT DATA,8X,5HNDP =,I3///
1 30H I XD YD ZD /)
2021 FORMAT(1X)
2022 FORMAT(5X,I2,2X,3F7.2)
2030 FORMAT(1H1,6HTTIDBS/////3X,17HIDBVIP SUBROUTINE///
1 26X,10HZI1(XI,YI))
2031 FORMAT(7X,2HXI,4X,3HYI=/12X,5F7.2/)
2032 FORMAT(1X/1X,F9.2,2X,5F7.2)
2040 FORMAT(1X/////3X,10HDIFFERENCE///
1 25X,11HDZI1(XI,YI))
2050 FORMAT(1H1,6HTTIDBS/////3X,17HIDSFFT SUBROUTINE///
1 26X,10HZI2(XI,YI))
2060 FORMAT(1X/////3X,10HDIFFERENCE///
1 25X,11HDZI2(XI,YI))
END
SUBROUTINE IDBVIP(MD,NCP,NDP,XD,YD,ZD,NIP,XI,YI,ZI,
1 IWK,WK)
C THIS SUBROUTINE PERFORMS BIVARIATE INTERPOLATION WHEN THE PRO-
C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
C DISTRIBUTED IN THE PLANE.
C THE INPUT PARAMETERS ARE
C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
C = 1 FOR NEW NCP AND/OR NEW XD-YD,
C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT
C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
C XD = ARRAY OF DIMENSION NDP CONTAINING THE X
C COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y
C COORDINATES OF THE DATA POINTS,
C ZD = ARRAY OF DIMENSION NDP CONTAINING THE Z
C COORDINATES OF THE DATA POINTS,
C NIP = NUMBER OF OUTPUT POINTS AT WHICH INTERPOLATION
C IS TO BE PERFORMED (MUST BE 1 OR GREATER),
C XI = ARRAY OF DIMENSION NIP CONTAINING THE X
C COORDINATES OF THE OUTPUT POINTS,
C YI = ARRAY OF DIMENSION NIP CONTAINING THE Y
C COORDINATES OF THE OUTPUT POINTS.
C THE OUTPUT PARAMETER IS
C ZI = ARRAY OF DIMENSION NIP WHERE INTERPOLATED Z
C VALUES ARE TO BE STORED.
C THE OTHER PARAMETERS ARE
C IWK = INTEGER ARRAY OF DIMENSION
C MAX0(31,27+NCP)*NDP+NIP
C USED INTERNALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A
C WORK AREA.
C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
C YD ARRAYS MUST BE MADE WITH MD=1. THE CALL WITH MD=2 MUST BE
C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS. THE CALL WITH
C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
C AND NIP VALUES AND WITH THE SAME CONTENTS OF THE XD, YD, XI,
C AND YI ARRAYS. BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C THIS SUBROUTINE CALLS THE IDCLDP, IDLCTN, IDPDRV, IDPTIP, AND
C IDTANG SUBROUTINES.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),XI(1000),YI(1000),
1 ZI(1000),IWK(4100),WK(800)
COMMON/IDLC/NIT
COMMON/IDPI/ITPV
DATA LUN/6/
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
C (FOR MD=1,2,3)
10 MD0=MD
NCP0=NCP
NDP0=NDP
NIP0=NIP
C ERROR CHECK. (FOR MD=1,2,3)
20 IF(MD0.LT.1.OR.MD0.GT.3) GO TO 90
IF(NCP0.LT.2.OR.NCP0.GE.NDP0) GO TO 90
IF(NDP0.LT.4) GO TO 90
IF(NIP0.LT.1) GO TO 90
IF(MD0.GE.2) GO TO 21
IWK(1)=NCP0
IWK(2)=NDP0
GO TO 22
21 NCPPV=IWK(1)
NDPPV=IWK(2)
IF(NCP0.NE.NCPPV) GO TO 90
IF(NDP0.NE.NDPPV) GO TO 90
22 IF(MD0.GE.3) GO TO 23
IWK(3)=NIP
GO TO 30
23 NIPPV=IWK(3)
IF(NIP0.NE.NIPPV) GO TO 90
C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY. (FOR MD=1,2,3)
30 JWIPT=16
JWIWL=6*NDP0+1
JWIWK=JWIWL
JWIPL=24*NDP0+1
JWIWP=30*NDP0+1
JWIPC=27*NDP0+1
JWIT0=MAX0(31,27+NCP0)*NDP0
C TRIANGULATES THE X-Y PLANE. (FOR MD=1)
40 IF(MD0.GT.1) GO TO 50
CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
1 IWK(JWIWL),IWK(JWIWP),WK)
IWK(5)=NT
IWK(6)=NL
IF(NT.EQ.0) RETURN
C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1)
50 IF(MD0.GT.1) GO TO 60
CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC))
IF(IWK(JWIPC).EQ.0) RETURN
C LOCATES ALL POINTS AT WHICH INTERPOLATION IS TO BE PERFORMED.
C (FOR MD=1,2)
60 IF(MD0.EQ.3) GO TO 70
NIT=0
JWIT=JWIT0
DO 61 IIP=1,NIP0
JWIT=JWIT+1
CALL IDLCTN(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
1 XI(IIP),YI(IIP),IWK(JWIT),IWK(JWIWK),WK)
61 CONTINUE
C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
C (FOR MD=1,2,3)
70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3)
80 ITPV=0
JWIT=JWIT0
DO 81 IIP=1,NIP0
JWIT=JWIT+1
CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
1 IWK(JWIT),XI(IIP),YI(IIP),ZI(IIP))
81 CONTINUE
RETURN
C ERROR EXIT
90 WRITE (LUN,2090) MD0,NCP0,NDP0,NIP0
RETURN
C FORMAT STATEMENT FOR ERROR MESSAGE
2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S)./
1 7H MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6,
2 10X,5HNIP =,I6/
3 35H ERROR DETECTED IN ROUTINE IDBVIP/)
END
SUBROUTINE IDCLDP(NDP,XD,YD,NCP,IPC)
C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST
C TO EACH OF THE DATA POINT.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C COORDINATES OF THE DATA POINTS,
C NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA
C POINTS.
C THE OUTPUT PARAMETER IS
C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE
C POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C EACH OF THE NDP DATA POINTS ARE TO BE STORED.
C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST
C NOT EXCEED 25.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),IPC(400)
DIMENSION DSQ0(25),IPC0(25)
DATA NCPMX/25/, LUN/6/
C STATEMENT FUNCTION
DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
C PRELIMINARY PROCESSING
10 NDP0=NDP
NCP0=NCP
IF(NDP0.LT.2) GO TO 90
IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0) GO TO 90
C CALCULATION
20 DO 59 IP1=1,NDP0
C - SELECTS NCP POINTS.
X1=XD(IP1)
Y1=YD(IP1)
J1=0
DSQMX=0.0
DO 22 IP2=1,NDP0
IF(IP2.EQ.IP1) GO TO 22
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
J1=J1+1
DSQ0(J1)=DSQI
IPC0(J1)=IP2
IF(DSQI.LE.DSQMX) GO TO 21
DSQMX=DSQI
JMX=J1
21 IF(J1.GE.NCP0) GO TO 23
22 CONTINUE
23 IP2MN=IP2+1
IF(IP2MN.GT.NDP0) GO TO 30
DO 25 IP2=IP2MN,NDP0
IF(IP2.EQ.IP1) GO TO 25
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
IF(DSQI.GE.DSQMX) GO TO 25
DSQ0(JMX)=DSQI
IPC0(JMX)=IP2
DSQMX=0.0
DO 24 J1=1,NCP0
IF(DSQ0(J1).LE.DSQMX) GO TO 24
DSQMX=DSQ0(J1)
JMX=J1
24 CONTINUE
25 CONTINUE
C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR.
30 IP2=IPC0(1)
DX12=XD(IP2)-X1
DY12=YD(IP2)-Y1
DO 31 J3=2,NCP0
IP3=IPC0(J3)
DX13=XD(IP3)-X1
DY13=YD(IP3)-Y1
IF((DY13*DX12-DX13*DY12).NE.0.0) GO TO 50
31 CONTINUE
C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT.
40 NCLPT=0
DO 43 IP3=1,NDP0
IF(IP3.EQ.IP1) GO TO 43
DO 41 J4=1,NCP0
IF(IP3.EQ.IPC0(J4)) GO TO 43
41 CONTINUE
DX13=XD(IP3)-X1
DY13=YD(IP3)-Y1
IF((DY13*DX12-DX13*DY12).EQ.0.0) GO TO 43
DSQI=DSQF(X1,Y1,XD(IP3),YD(IP3))
IF(NCLPT.EQ.0) GO TO 42
IF(DSQI.GE.DSQMN) GO TO 43
42 NCLPT=1
DSQMN=DSQI
IP3MN=IP3
43 CONTINUE
IF(NCLPT.EQ.0) GO TO 91
DSQMX=DSQMN
IPC0(JMX)=IP3MN
C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY.
50 J1=(IP1-1)*NCP0
DO 51 J2=1,NCP0
J1=J1+1
IPC(J1)=IPC0(J2)
51 CONTINUE
59 CONTINUE
RETURN
C ERROR EXIT
90 WRITE (LUN,2090)
GO TO 92
91 WRITE (LUN,2091)
92 WRITE (LUN,2092) NDP0,NCP0
IPC(1)=0
RETURN
C FORMAT STATEMENTS FOR ERROR MESSAGES
2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S).)
2091 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS.)
2092 FORMAT(8H NDP =,I5,5X,5HNCP =,I5/
1 35H ERROR DETECTED IN ROUTINE IDCLDP/)
END
SUBROUTINE IDGRID(XD, YD, NT, IPT, NL, IPL, NXI, NYI, XI, YI,
* NGP, IGP)
C THIS SUBROUTINE ORGANIZES GRID POINTS FOR SURFACE FITTING BY
C SORTING THEM IN ASCENDING ORDER OF TRIANGLE NUMBERS AND OF THE
C BORDER LINE SEGMENT NUMBER.
C THE INPUT PARAMETERS ARE
C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C COORDINATES OF THE DATA POINTS, WHERE NDP IS THE
C NUMBER OF THE DATA POINTS,
C NT = NUMBER OF TRIANGLES,
C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
C POINT NUMBERS OF THE END POINTS OF THE BORDER
C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
C NUMBERS,
C NXI = NUMBER OF GRID POINTS IN THE X COORDINATE,
C NYI = NUMBER OF GRID POINTS IN THE Y COORDINATE,
C XI,YI = ARRAYS OF DIMENSION NXI AND NYI CONTAINING
C THE X AND Y COORDINATES OF THE GRID POINTS,
C RESPECTIVELY.
C THE OUTPUT PARAMETERS ARE
C NGP = INTEGER ARRAY OF DIMENSION 2*(NT+2*NL) WHERE THE
C NUMBER OF GRID POINTS THAT BELONG TO EACH OF THE
C TRIANGLES OR OF THE BORDER LINE SEGMENTS ARE TO
C BE STORED,
C IGP = INTEGER ARRAY OF DIMENSION NXI*NYI WHERE THE
C GRID POINT NUMBERS ARE TO BE STORED IN ASCENDING
C ORDER OF THE TRIANGLE NUMBER AND THE BORDER LINE
C SEGMENT NUMBER.
C DECLARATION STATEMENTS
DIMENSION XD(100), YD(100), IPT(585), IPL(300), XI(101),
* YI(101), NGP(800), IGP(10201)
C STATEMENT FUNCTIONS
SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3)
SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2)
C PRELIMINARY PROCESSING
NT0 = NT
NL0 = NL
NXI0 = NXI
NYI0 = NYI
NXINYI = NXI0*NYI0
XIMN = AMIN1(XI(1),XI(NXI0))
XIMX = AMAX1(XI(1),XI(NXI0))
YIMN = AMIN1(YI(1),YI(NYI0))
YIMX = AMAX1(YI(1),YI(NYI0))
C DETERMINES GRID POINTS INSIDE THE DATA AREA.
JNGP0 = 0
JNGP1 = 2*(NT0+2*NL0) + 1
JIGP0 = 0
JIGP1 = NXINYI + 1
DO 160 IT0=1,NT0
NGP0 = 0
NGP1 = 0
IT0T3 = IT0*3
IP1 = IPT(IT0T3-2)
IP2 = IPT(IT0T3-1)
IP3 = IPT(IT0T3)
X1 = XD(IP1)
Y1 = YD(IP1)
X2 = XD(IP2)
Y2 = YD(IP2)
X3 = XD(IP3)
Y3 = YD(IP3)
XMN = AMIN1(X1,X2,X3)
XMX = AMAX1(X1,X2,X3)
YMN = AMIN1(Y1,Y2,Y3)
YMX = AMAX1(Y1,Y2,Y3)
INSD = 0
DO 20 IXI=1,NXI0
IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 10
IF (INSD.EQ.0) GO TO 20
IXIMX = IXI - 1
GO TO 30
10 IF (INSD.EQ.1) GO TO 20
INSD = 1
IXIMN = IXI
20 CONTINUE
IF (INSD.EQ.0) GO TO 150
IXIMX = NXI0
30 DO 140 IYI=1,NYI0
YII = YI(IYI)
IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 140
DO 130 IXI=IXIMN,IXIMX
XII = XI(IXI)
L = 0
IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 130, 40, 50
40 L = 1
50 IF (SIDE(X2,Y2,X3,Y3,XII,YII)) 130, 60, 70
60 L = 1
70 IF (SIDE(X3,Y3,X1,Y1,XII,YII)) 130, 80, 90
80 L = 1
90 IZI = NXI0*(IYI-1) + IXI
IF (L.EQ.1) GO TO 100
NGP0 = NGP0 + 1
JIGP0 = JIGP0 + 1
IGP(JIGP0) = IZI
GO TO 130
100 IF (JIGP1.GT.NXINYI) GO TO 120
DO 110 JIGP1I=JIGP1,NXINYI
IF (IZI.EQ.IGP(JIGP1I)) GO TO 130
110 CONTINUE
120 NGP1 = NGP1 + 1
JIGP1 = JIGP1 - 1
IGP(JIGP1) = IZI
130 CONTINUE
140 CONTINUE
150 JNGP0 = JNGP0 + 1
NGP(JNGP0) = NGP0
JNGP1 = JNGP1 - 1
NGP(JNGP1) = NGP1
160 CONTINUE
C DETERMINES GRID POINTS OUTSIDE THE DATA AREA.
C - IN SEMI-INFINITE RECTANGULAR AREA.
DO 450 IL0=1,NL0
NGP0 = 0
NGP1 = 0
IL0T3 = IL0*3
IP1 = IPL(IL0T3-2)
IP2 = IPL(IL0T3-1)
X1 = XD(IP1)
Y1 = YD(IP1)
X2 = XD(IP2)
Y2 = YD(IP2)
XMN = XIMN
XMX = XIMX
YMN = YIMN
YMX = YIMX
IF (Y2.GE.Y1) XMN = AMIN1(X1,X2)
IF (Y2.LE.Y1) XMX = AMAX1(X1,X2)
IF (X2.LE.X1) YMN = AMIN1(Y1,Y2)
IF (X2.GE.X1) YMX = AMAX1(Y1,Y2)
INSD = 0
DO 180 IXI=1,NXI0
IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 170
IF (INSD.EQ.0) GO TO 180
IXIMX = IXI - 1
GO TO 190
170 IF (INSD.EQ.1) GO TO 180
INSD = 1
IXIMN = IXI
180 CONTINUE
IF (INSD.EQ.0) GO TO 310
IXIMX = NXI0
190 DO 300 IYI=1,NYI0
YII = YI(IYI)
IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 300
DO 290 IXI=IXIMN,IXIMX
XII = XI(IXI)
L = 0
IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 210, 200, 290
200 L = 1
210 IF (SPDT(X2,Y2,X1,Y1,XII,YII)) 290, 220, 230
220 L = 1
230 IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 290, 240, 250
240 L = 1
250 IZI = NXI0*(IYI-1) + IXI
IF (L.EQ.1) GO TO 260
NGP0 = NGP0 + 1
JIGP0 = JIGP0 + 1
IGP(JIGP0) = IZI
GO TO 290
260 IF (JIGP1.GT.NXINYI) GO TO 280
DO 270 JIGP1I=JIGP1,NXINYI
IF (IZI.EQ.IGP(JIGP1I)) GO TO 290
270 CONTINUE
280 NGP1 = NGP1 + 1
JIGP1 = JIGP1 - 1
IGP(JIGP1) = IZI
290 CONTINUE
300 CONTINUE
310 JNGP0 = JNGP0 + 1
NGP(JNGP0) = NGP0
JNGP1 = JNGP1 - 1
NGP(JNGP1) = NGP1
C - IN SEMI-INFINITE TRIANGULAR AREA.
NGP0 = 0
NGP1 = 0
ILP1 = MOD(IL0,NL0) + 1
ILP1T3 = ILP1*3
IP3 = IPL(ILP1T3-1)
X3 = XD(IP3)
Y3 = YD(IP3)
XMN = XIMN
XMX = XIMX
YMN = YIMN
YMX = YIMX
IF (Y3.GE.Y2 .AND. Y2.GE.Y1) XMN = X2
IF (Y3.LE.Y2 .AND. Y2.LE.Y1) XMX = X2
IF (X3.LE.X2 .AND. X2.LE.X1) YMN = Y2
IF (X3.GE.X2 .AND. X2.GE.X1) YMX = Y2
INSD = 0
DO 330 IXI=1,NXI0
IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 320
IF (INSD.EQ.0) GO TO 330
IXIMX = IXI - 1
GO TO 340
320 IF (INSD.EQ.1) GO TO 330
INSD = 1
IXIMN = IXI
330 CONTINUE
IF (INSD.EQ.0) GO TO 440
IXIMX = NXI0
340 DO 430 IYI=1,NYI0
YII = YI(IYI)
IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 430
DO 420 IXI=IXIMN,IXIMX
XII = XI(IXI)
L = 0
IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 360, 350, 420
350 L = 1
360 IF (SPDT(X3,Y3,X2,Y2,XII,YII)) 380, 370, 420
370 L = 1
380 IZI = NXI0*(IYI-1) + IXI
IF (L.EQ.1) GO TO 390
NGP0 = NGP0 + 1
JIGP0 = JIGP0 + 1
IGP(JIGP0) = IZI
GO TO 420
390 IF (JIGP1.GT.NXINYI) GO TO 410
DO 400 JIGP1I=JIGP1,NXINYI
IF (IZI.EQ.IGP(JIGP1I)) GO TO 420
400 CONTINUE
410 NGP1 = NGP1 + 1
JIGP1 = JIGP1 - 1
IGP(JIGP1) = IZI
420 CONTINUE
430 CONTINUE
440 JNGP0 = JNGP0 + 1
NGP(JNGP0) = NGP0
JNGP1 = JNGP1 - 1
NGP(JNGP1) = NGP1
450 CONTINUE
RETURN
END
SUBROUTINE IDLCTN(NDP, XD, YD, NT, IPT, NL, IPL, XII, YII, ITI,
* IWK, WK)
C THIS SUBROUTINE LOCATES A POINT, I.E., DETERMINES TO WHAT TRI-
C ANGLE A GIVEN POINT (XII,YII) BELONGS. WHEN THE GIVEN POINT
C DOES NOT LIE INSIDE THE DATA AREA, THIS SUBROUTINE DETERMINES
C THE BORDER LINE SEGMENT WHEN THE POINT LIES IN AN OUTSIDE
C RECTANGULAR AREA, AND TWO BORDER LINE SEGMENTS WHEN THE POINT
C LIES IN AN OUTSIDE TRIANGULAR AREA.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C COORDINATES OF THE DATA POINTS,
C NT = NUMBER OF TRIANGLES,
C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
C POINT NUMBERS OF THE END POINTS OF THE BORDER
C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
C NUMBERS,
C XII,YII = X AND Y COORDINATES OF THE POINT TO BE
C LOCATED.
C THE OUTPUT PARAMETER IS
C ITI = TRIANGLE NUMBER, WHEN THE POINT IS INSIDE THE
C DATA AREA, OR
C TWO BORDER LINE SEGMENT NUMBERS, IL1 AND IL2,
C CODED TO IL1*(NT+NL)+IL2, WHEN THE POINT IS
C OUTSIDE THE DATA AREA.
C THE OTHER PARAMETERS ARE
C IWK = INTEGER ARRAY OF DIMENSION 18*NDP USED INTER-
C NALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A
C WORK AREA.
C DECLARATION STATEMENTS
DIMENSION XD(100), YD(100), IPT(585), IPL(300), IWK(1800),
* WK(800)
DIMENSION NTSC(9), IDSC(9)
COMMON /IDLC/ NIT
C STATEMENT FUNCTIONS
SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3)
SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2)
C PRELIMINARY PROCESSING
NDP0 = NDP
NT0 = NT
NL0 = NL
NTL = NT0 + NL0
X0 = XII
Y0 = YII
C PROCESSING FOR A NEW SET OF DATA POINTS
IF (NIT.NE.0) GO TO 80
NIT = 1
C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS.
XMN = XD(1)
XMX = XMN
YMN = YD(1)
YMX = YMN
DO 10 IDP=2,NDP0
XI = XD(IDP)
YI = YD(IDP)
XMN = AMIN1(XI,XMN)
XMX = AMAX1(XI,XMX)
YMN = AMIN1(YI,YMN)
YMX = AMAX1(YI,YMX)
10 CONTINUE
XS1 = (XMN+XMN+XMX)/3.0
XS2 = (XMN+XMX+XMX)/3.0
YS1 = (YMN+YMN+YMX)/3.0
YS2 = (YMN+YMX+YMX)/3.0
C - DETERMINES AND STORES IN THE IWK ARRAY TRIANGLE NUMBERS OF
C - THE TRIANGLES ASSOCIATED WITH EACH OF THE NINE SECTIONS.
DO 20 ISC=1,9
NTSC(ISC) = 0
IDSC(ISC) = 0
20 CONTINUE
IT0T3 = 0
JWK = 0
DO 70 IT0=1,NT0
IT0T3 = IT0T3 + 3
I1 = IPT(IT0T3-2)
I2 = IPT(IT0T3-1)
I3 = IPT(IT0T3)
XMN = AMIN1(XD(I1),XD(I2),XD(I3))
XMX = AMAX1(XD(I1),XD(I2),XD(I3))
YMN = AMIN1(YD(I1),YD(I2),YD(I3))
YMX = AMAX1(YD(I1),YD(I2),YD(I3))
IF (YMN.GT.YS1) GO TO 30
IF (XMN.LE.XS1) IDSC(1) = 1
IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(2) = 1
IF (XMX.GE.XS2) IDSC(3) = 1
30 IF (YMX.LT.YS1 .OR. YMN.GT.YS2) GO TO 40
IF (XMN.LE.XS1) IDSC(4) = 1
IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(5) = 1
IF (XMX.GE.XS2) IDSC(6) = 1
40 IF (YMX.LT.YS2) GO TO 50
IF (XMN.LE.XS1) IDSC(7) = 1
IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(8) = 1
IF (XMX.GE.XS2) IDSC(9) = 1
50 DO 60 ISC=1,9
IF (IDSC(ISC).EQ.0) GO TO 60
JIWK = 9*NTSC(ISC) + ISC
IWK(JIWK) = IT0
NTSC(ISC) = NTSC(ISC) + 1
IDSC(ISC) = 0
60 CONTINUE
C - STORES IN THE WK ARRAY THE MINIMUM AND MAXIMUM OF THE X AND
C - Y COORDINATE VALUES FOR EACH OF THE TRIANGLE.
JWK = JWK + 4
WK(JWK-3) = XMN
WK(JWK-2) = XMX
WK(JWK-1) = YMN
WK(JWK) = YMX
70 CONTINUE
GO TO 110
C CHECKS IF IN THE SAME TRIANGLE AS PREVIOUS.
80 IT0 = ITIPV
IF (IT0.GT.NT0) GO TO 90
IT0T3 = IT0*3
IP1 = IPT(IT0T3-2)
X1 = XD(IP1)
Y1 = YD(IP1)
IP2 = IPT(IT0T3-1)
X2 = XD(IP2)
Y2 = YD(IP2)
IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110
IP3 = IPT(IT0T3)
X3 = XD(IP3)
Y3 = YD(IP3)
IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 110
IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 110
GO TO 170
C CHECKS IF ON THE SAME BORDER LINE SEGMENT.
90 IL1 = IT0/NTL
IL2 = IT0 - IL1*NTL
IL1T3 = IL1*3
IP1 = IPL(IL1T3-2)
X1 = XD(IP1)
Y1 = YD(IP1)
IP2 = IPL(IL1T3-1)
X2 = XD(IP2)
Y2 = YD(IP2)
IF (IL2.NE.IL1) GO TO 100
IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110
IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 110
IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110
GO TO 170
C CHECKS IF BETWEEN THE SAME TWO BORDER LINE SEGMENTS.
100 IF (SPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110
IP3 = IPL(3*IL2-1)
X3 = XD(IP3)
Y3 = YD(IP3)
IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 170
C LOCATES INSIDE THE DATA AREA.
C - DETERMINES THE SECTION IN WHICH THE POINT IN QUESTION LIES.
110 ISC = 1
IF (X0.GE.XS1) ISC = ISC + 1
IF (X0.GE.XS2) ISC = ISC + 1
IF (Y0.GE.YS1) ISC = ISC + 3
IF (Y0.GE.YS2) ISC = ISC + 3
C - SEARCHES THROUGH THE TRIANGLES ASSOCIATED WITH THE SECTION.
NTSCI = NTSC(ISC)
IF (NTSCI.LE.0) GO TO 130
JIWK = -9 + ISC
DO 120 ITSC=1,NTSCI
JIWK = JIWK + 9
IT0 = IWK(JIWK)
JWK = IT0*4
IF (X0.LT.WK(JWK-3)) GO TO 120
IF (X0.GT.WK(JWK-2)) GO TO 120
IF (Y0.LT.WK(JWK-1)) GO TO 120
IF (Y0.GT.WK(JWK)) GO TO 120
IT0T3 = IT0*3
IP1 = IPT(IT0T3-2)
X1 = XD(IP1)
Y1 = YD(IP1)
IP2 = IPT(IT0T3-1)
X2 = XD(IP2)
Y2 = YD(IP2)
IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 120
IP3 = IPT(IT0T3)
X3 = XD(IP3)
Y3 = YD(IP3)
IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 120
IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 120
GO TO 170
120 CONTINUE
C LOCATES OUTSIDE THE DATA AREA.
130 DO 150 IL1=1,NL0
IL1T3 = IL1*3
IP1 = IPL(IL1T3-2)
X1 = XD(IP1)
Y1 = YD(IP1)
IP2 = IPL(IL1T3-1)
X2 = XD(IP2)
Y2 = YD(IP2)
IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 150
IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 140
IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 150
IL2 = IL1
GO TO 160
140 IL2 = MOD(IL1,NL0) + 1
IP3 = IPL(3*IL2-1)
X3 = XD(IP3)
Y3 = YD(IP3)
IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 160
150 CONTINUE
IT0 = 1
GO TO 170
160 IT0 = IL1*NTL + IL2
C NORMAL EXIT
170 ITI = IT0
ITIPV = IT0
RETURN
END
SUBROUTINE IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD)
C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND
C SECOND ORDER AT THE DATA POINTS.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
C Y, AND Z COORDINATES OF THE DATA POINTS,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT,
C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING
C THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C EACH OF THE NDP DATA POINTS.
C THE OUTPUT PARAMETER IS
C PD = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED
C ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA
C POINTS ARE TO BE STORED.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),IPC(400),PD(500)
REAL NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY
C PRELIMINARY PROCESSING
10 NDP0=NDP
NCP0=NCP
NCPM1=NCP0-1
C ESTIMATION OF ZX AND ZY
20 DO 24 IP0=1,NDP0
X0=XD(IP0)
Y0=YD(IP0)
Z0=ZD(IP0)
NMX=0.0
NMY=0.0
NMZ=0.0
JIPC0=NCP0*(IP0-1)
DO 23 IC1=1,NCPM1
JIPC=JIPC0+IC1
IPI=IPC(JIPC)
DX1=XD(IPI)-X0
DY1=YD(IPI)-Y0
DZ1=ZD(IPI)-Z0
IC2MN=IC1+1
DO 22 IC2=IC2MN,NCP0
JIPC=JIPC0+IC2
IPI=IPC(JIPC)
DX2=XD(IPI)-X0
DY2=YD(IPI)-Y0
DNMZ=DX1*DY2-DY1*DX2
IF(DNMZ.EQ.0.0) GO TO 22
DZ2=ZD(IPI)-Z0
DNMX=DY1*DZ2-DZ1*DY2
DNMY=DZ1*DX2-DX1*DZ2
IF(DNMZ.GE.0.0) GO TO 21
DNMX=-DNMX
DNMY=-DNMY
DNMZ=-DNMZ
21 NMX=NMX+DNMX
NMY=NMY+DNMY
NMZ=NMZ+DNMZ
22 CONTINUE
23 CONTINUE
JPD0=5*IP0
PD(JPD0-4)=-NMX/NMZ
PD(JPD0-3)=-NMY/NMZ
24 CONTINUE
C ESTIMATION OF ZXX, ZXY, AND ZYY
30 DO 34 IP0=1,NDP0
JPD0=JPD0+5
X0=XD(IP0)
JPD0=5*IP0
Y0=YD(IP0)
ZX0=PD(JPD0-4)
ZY0=PD(JPD0-3)
NMXX=0.0
NMXY=0.0
NMYX=0.0
NMYY=0.0
NMZ =0.0
JIPC0=NCP0*(IP0-1)
DO 33 IC1=1,NCPM1
JIPC=JIPC0+IC1
IPI=IPC(JIPC)
DX1=XD(IPI)-X0
DY1=YD(IPI)-Y0
JPD=5*IPI
DZX1=PD(JPD-4)-ZX0
DZY1=PD(JPD-3)-ZY0
IC2MN=IC1+1
DO 32 IC2=IC2MN,NCP0
JIPC=JIPC0+IC2
IPI=IPC(JIPC)
DX2=XD(IPI)-X0
DY2=YD(IPI)-Y0
DNMZ =DX1*DY2 -DY1*DX2
IF(DNMZ.EQ.0.0) GO TO 32
JPD=5*IPI
DZX2=PD(JPD-4)-ZX0
DZY2=PD(JPD-3)-ZY0
DNMXX=DY1*DZX2-DZX1*DY2
DNMXY=DZX1*DX2-DX1*DZX2
DNMYX=DY1*DZY2-DZY1*DY2
DNMYY=DZY1*DX2-DX1*DZY2
IF(DNMZ.GE.0.0) GO TO 31
DNMXX=-DNMXX
DNMXY=-DNMXY
DNMYX=-DNMYX
DNMYY=-DNMYY
DNMZ =-DNMZ
31 NMXX=NMXX+DNMXX
NMXY=NMXY+DNMXY
NMYX=NMYX+DNMYX
ID010010
NMYY=NMYY+DNMYY
NMZ =NMZ +DNMZ
32 CONTINUE
33 CONTINUE
PD(JPD0-2)=-NMXX/NMZ
PD(JPD0-1)=-(NMXY+NMYX)/(2.0*NMZ)
PD(JPD0) =-NMYY/NMZ
34 CONTINUE
RETURN
END
SUBROUTINE IDPTIP(XD,YD,ZD,NT,IPT,NL,IPL,PDD,ITI,XII,YII,
1 ZII)
C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPOLA-
C TION, I.E., DETERMINES THE Z VALUE AT A POINT.
C THE INPUT PARAMETERS ARE
C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
C Y, AND Z COORDINATES OF THE DATA POINTS, WHERE
C NDP IS THE NUMBER OF THE DATA POINTS,
C NT = NUMBER OF TRIANGLES,
C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
C POINT NUMBERS OF THE END POINTS OF THE BORDER
C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
C NUMBERS,
C PDD = ARRAY OF DIMENSION 5*NDP CONTAINING THE PARTIAL
C DERIVATIVES AT THE DATA POINTS,
C ITI = TRIANGLE NUMBER OF THE TRIANGLE IN WHICH LIES
C THE POINT FOR WHICH INTERPOLATION IS TO BE
C PERFORMED,
C XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH
C INTERPOLATION IS TO BE PERFORMED.
C THE OUTPUT PARAMETER IS
C ZII = INTERPOLATED Z VALUE.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),IPT(585),IPL(300),
1 PDD(500)
COMMON/IDPI/ITPV
DIMENSION X(3),Y(3),Z(3),PD(15),
1 ZU(3),ZV(3),ZUU(3),ZUV(3),ZVV(3)
REAL LU,LV
EQUIVALENCE (P5,P50)
C PRELIMINARY PROCESSING
10 IT0=ITI
NTL=NT+NL
IF(IT0.LE.NTL) GO TO 20
IL1=IT0/NTL
IL2=IT0-IL1*NTL
IF(IL1.EQ.IL2) GO TO 40
GO TO 60
C CALCULATION OF ZII BY INTERPOLATION.
C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
20 IF(IT0.EQ.ITPV) GO TO 30
C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE
C VERTEXES.
21 JIPT=3*(IT0-1)
JPD=0
DO 23 I=1,3
JIPT=JIPT+1
IDP=IPT(JIPT)
X(I)=XD(IDP)
Y(I)=YD(IDP)
Z(I)=ZD(IDP)
JPDD=5*(IDP-1)
DO 22 KPD=1,5
JPD=JPD+1
JPDD=JPDD+1
PD(JPD)=PDD(JPDD)
22 CONTINUE
23 CONTINUE
C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
C AND VICE VERSA.
24 X0=X(1)
Y0=Y(1)
A=X(2)-X0
B=X(3)-X0
C=Y(2)-Y0
D=Y(3)-Y0
AD=A*D
BC=B*C
DLT=AD-BC
AP= D/DLT
BP=-B/DLT
CP=-C/DLT
DP= A/DLT
C CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE
C TRIANGLE FOR THE U-V COORDINATE SYSTEM.
25 AA=A*A
ACT2=2.0*A*C
CC=C*C
AB=A*B
ADBC=AD+BC
CD=C*D
BB=B*B
BDT2=2.0*B*D
DD=D*D
DO 26 I=1,3
JPD=5*I
ZU(I)=A*PD(JPD-4)+C*PD(JPD-3)
ZV(I)=B*PD(JPD-4)+D*PD(JPD-3)
ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD)
ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD)
ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD)
26 CONTINUE
C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
27 P00=Z(1)
P10=ZU(1)
P01=ZV(1)
P20=0.5*ZUU(1)
P11=ZUV(1)
P02=0.5*ZVV(1)
H1=Z(2)-P00-P10-P20
H2=ZU(2)-P10-ZUU(1)
H3=ZUU(2)-ZUU(1)
P30= 10.0*H1-4.0*H2+0.5*H3
P40=-15.0*H1+7.0*H2 -H3
P50= 6.0*H1-3.0*H2+0.5*H3
H1=Z(3)-P00-P01-P02
H2=ZV(3)-P01-ZVV(1)
H3=ZVV(3)-ZVV(1)
P03= 10.0*H1-4.0*H2+0.5*H3
P04=-15.0*H1+7.0*H2 -H3
P05= 6.0*H1-3.0*H2+0.5*H3
LU=SQRT(AA+CC)
LV=SQRT(BB+DD)
THXU=ATAN2(C,A)
THUV=ATAN2(D,B)-THXU
CSUV=COS(THUV)
P41=5.0*LV*CSUV/LU*P50
P14=5.0*LU*CSUV/LV*P05
H1=ZV(2)-P01-P11-P41
H2=ZUV(2)-P11-4.0*P41
P21= 3.0*H1-H2
P31=-2.0*H1+H2
H1=ZU(3)-P10-P11-P14
H2=ZUV(3)-P11-4.0*P14
P12= 3.0*H1-H2
P13=-2.0*H1+H2
THUS=ATAN2(D-C,B-A)-THXU
THSV=THUV-THUS
AA= SIN(THSV)/LU
BB=-COS(THSV)/LU
CC= SIN(THUS)/LV
DD= COS(THUS)/LV
AC=AA*CC
AD=AA*DD
BC=BB*CC
G1=AA*AC*(3.0*BC+2.0*AD)
G2=CC*AC*(3.0*AD+2.0*BC)
H1=-AA*AA*AA*(5.0*AA*BB*P50+(4.0*BC+AD)*P41)
1 -CC*CC*CC*(5.0*CC*DD*P05+(4.0*AD+BC)*P14)
H2=0.5*ZVV(2)-P02-P12
H3=0.5*ZUU(3)-P20-P21
P22=(G1*H2+G2*H3-H1)/(G1+G2)
P32=H2-P22
P23=H3-P22
ITPV=IT0
C CONVERTS XII AND YII TO U-V SYSTEM.
30 DX=XII-X0
DY=YII-Y0
U=AP*DX+BP*DY
V=CP*DX+DP*DY
C EVALUATES THE POLYNOMIAL.
31 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05))))
P1=P10+V*(P11+V*(P12+V*(P13+V*P14)))
P2=P20+V*(P21+V*(P22+V*P23))
P3=P30+V*(P31+V*P32)
P4=P40+V*P41
ZII=P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5))))
RETURN
C CALCULATION OF ZII BY EXTRAPOLATION IN THE RECTANGLE.
C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
40 IF(IT0.EQ.ITPV) GO TO 50
C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE END
C POINTS OF THE BORDER LINE SEGMENT.
41 JIPL=3*(IL1-1)
JPD=0
DO 43 I=1,2
JIPL=JIPL+1
IDP=IPL(JIPL)
X(I)=XD(IDP)
Y(I)=YD(IDP)
Z(I)=ZD(IDP)
JPDD=5*(IDP-1)
DO 42 KPD=1,5
JPD=JPD+1
JPDD=JPDD+1
PD(JPD)=PDD(JPDD)
42 CONTINUE
43 CONTINUE
C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
C AND VICE VERSA.
44 X0=X(1)
Y0=Y(1)
A=Y(2)-Y(1)
B=X(2)-X(1)
C=-B
D=A
AD=A*D
BC=B*C
DLT=AD-BC
AP= D/DLT
BP=-B/DLT
CP=-BP
DP= AP
C CONVERTS THE PARTIAL DERIVATIVES AT THE END POINTS OF THE
C BORDER LINE SEGMENT FOR THE U-V COORDINATE SYSTEM.
45 AA=A*A
ACT2=2.0*A*C
CC=C*C
AB=A*B
ADBC=AD+BC
CD=C*D
BB=B*B
BDT2=2.0*B*D
DD=D*D
DO 46 I=1,2
JPD=5*I
ZU(I)=A*PD(JPD-4)+C*PD(JPD-3)
ZV(I)=B*PD(JPD-4)+D*PD(JPD-3)
ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD)
ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD)
ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD)
46 CONTINUE
C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
47 P00=Z(1)
P10=ZU(1)
P01=ZV(1)
P20=0.5*ZUU(1)
P11=ZUV(1)
P02=0.5*ZVV(1)
H1=Z(2)-P00-P01-P02
H2=ZV(2)-P01-ZVV(1)
H3=ZVV(2)-ZVV(1)
P03= 10.0*H1-4.0*H2+0.5*H3
P04=-15.0*H1+7.0*H2 -H3
P05= 6.0*H1-3.0*H2+0.5*H3
H1=ZU(2)-P10-P11
H2=ZUV(2)-P11
P12= 3.0*H1-H2
P13=-2.0*H1+H2
P21=0.0
P23=-ZUU(2)+ZUU(1)
P22=-1.5*P23
ITPV=IT0
C CONVERTS XII AND YII TO U-V SYSTEM.
50 DX=XII-X0
DY=YII-Y0
U=AP*DX+BP*DY
V=CP*DX+DP*DY
C EVALUATES THE POLYNOMIAL.
51 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05))))
P1=P10+V*(P11+V*(P12+V*P13))
P2=P20+V*(P21+V*(P22+V*P23))
ZII=P0+U*(P1+U*P2)
RETURN
C CALCULATION OF ZII BY EXTRAPOLATION IN THE TRIANGLE.
C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
60 IF(IT0.EQ.ITPV) GO TO 70
C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE VERTEX
C OF THE TRIANGLE.
61 JIPL=3*IL2-2
IDP=IPL(JIPL)
X(1)=XD(IDP)
Y(1)=YD(IDP)
Z(1)=ZD(IDP)
JPDD=5*(IDP-1)
DO 62 KPD=1,5
JPDD=JPDD+1
PD(KPD)=PDD(JPDD)
62 CONTINUE
C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
67 P00=Z(1)
P10=PD(1)
P01=PD(2)
P20=0.5*PD(3)
P11=PD(4)
P02=0.5*PD(5)
ITPV=IT0
C CONVERTS XII AND YII TO U-V SYSTEM.
70 U=XII-X(1)
V=YII-Y(1)
C EVALUATES THE POLYNOMIAL.
71 P0=P00+V*(P01+V*P02)
P1=P10+V*P11
ZII=P0+U*(P1+U*P20)
RETURN
END
SUBROUTINE IDSFFT(MD,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI,
1 IWK,WK)
C THIS SUBROUTINE PERFORMS SMOOTH SURFACE FITTING WHEN THE PRO-
C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
C DISTRIBUTED IN THE PLANE.
C THE INPUT PARAMETERS ARE
C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
C = 1 FOR NEW NCP AND/OR NEW XD-YD,
C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT
C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
C XD = ARRAY OF DIMENSION NDP CONTAINING THE X
C COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y
C COORDINATES OF THE DATA POINTS,
C ZD = ARRAY OF DIMENSION NDP CONTAINING THE Z
C COORDINATES OF THE DATA POINTS,
C NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE
C (MUST BE 1 OR GREATER),
C NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE
C (MUST BE 1 OR GREATER),
C XI = ARRAY OF DIMENSION NXI CONTAINING THE X
C COORDINATES OF THE OUTPUT GRID POINTS,
C YI = ARRAY OF DIMENSION NYI CONTAINING THE Y
C COORDINATES OF THE OUTPUT GRID POINTS.
C THE OUTPUT PARAMETER IS
C ZI = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (NXI,NYI),
C WHERE THE INTERPOLATED Z VALUES AT THE OUTPUT
C GRID POINTS ARE TO BE STORED.
C THE OTHER PARAMETERS ARE
C IWK = INTEGER ARRAY OF DIMENSION
C MAX0(31,27+NCP)*NDP+NXI*NYI
C USED INTERNALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION 5*NDP USED INTERNALLY AS A
C WORK AREA.
C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
C YD ARRAYS MUST BE MADE WITH MD=1. THE CALL WITH MD=2 MUST BE
C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS. THE CALL WITH
C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
C NXI, AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD, YD,
C XI, AND YI ARRAYS. BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C THIS SUBROUTINE CALLS THE IDCLDP, IDGRID, IDPDRV, IDPTIP, AND
C IDTANG SUBROUTINES.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),XI(101),YI(101),
1 ZI(10201),IWK(13301),WK(500)
COMMON/IDPI/ITPV
DATA LUN/6/
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
C (FOR MD=1,2,3)
10 MD0=MD
NCP0=NCP
NDP0=NDP
NXI0=NXI
NYI0=NYI
C ERROR CHECK. (FOR MD=1,2,3)
20 IF(MD0.LT.1.OR.MD0.GT.3) GO TO 90
IF(NCP0.LT.2.OR.NCP0.GE.NDP0) GO TO 90
IF(NDP0.LT.4) GO TO 90
IF(NXI0.LT.1.OR.NYI0.LT.1) GO TO 90
IF(MD0.GE.2) GO TO 21
IWK(1)=NCP0
IWK(2)=NDP0
GO TO 22
21 NCPPV=IWK(1)
NDPPV=IWK(2)
IF(NCP0.NE.NCPPV) GO TO 90
IF(NDP0.NE.NDPPV) GO TO 90
22 IF(MD0.GE.3) GO TO 23
IWK(3)=NXI0
IWK(4)=NYI0
GO TO 30
23 NXIPV=IWK(3)
NYIPV=IWK(4)
IF(NXI0.NE.NXIPV) GO TO 90
IF(NYI0.NE.NYIPV) GO TO 90
C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY. (FOR MD=1,2,3)
30 JWIPT=16
JWIWL=6*NDP0+1
JWNGP0=JWIWL-1
JWIPL=24*NDP0+1
JWIWP=30*NDP0+1
JWIPC=27*NDP0+1
JWIGP0=MAX0(31,27+NCP0)*NDP0
C TRIANGULATES THE X-Y PLANE. (FOR MD=1)
40 IF(MD0.GT.1) GO TO 50
CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
1 IWK(JWIWL),IWK(JWIWP),WK)
IWK(5)=NT
IWK(6)=NL
IF(NT.EQ.0) RETURN
C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1)
50 IF(MD0.GT.1) GO TO 60
CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC))
IF(IWK(JWIPC).EQ.0) RETURN
C SORTS OUTPUT GRID POINTS IN ASCENDING ORDER OF THE TRIANGLE
C NUMBER AND THE BORDER LINE SEGMENT NUMBER. (FOR MD=1,2)
60 IF(MD0.EQ.3) GO TO 70
CALL IDGRID(XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),NXI0,NYI0,
1 XI,YI,IWK(JWNGP0+1),IWK(JWIGP0+1))
C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
C (FOR MD=1,2,3)
70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3)
80 ITPV=0
JIG0MX=0
JIG1MN=NXI0*NYI0+1
NNGP=NT+2*NL
DO 89 JNGP=1,NNGP
ITI=JNGP
IF(JNGP.LE.NT) GO TO 81
IL1=(JNGP-NT+1)/2
IL2=(JNGP-NT+2)/2
IF(IL2.GT.NL) IL2=1
ITI=IL1*(NT+NL)+IL2
81 JWNGP=JWNGP0+JNGP
NGP0=IWK(JWNGP)
IF(NGP0.EQ.0) GO TO 86
JIG0MN=JIG0MX+1
JIG0MX=JIG0MX+NGP0
DO 82 JIGP=JIG0MN,JIG0MX
JWIGP=JWIGP0+JIGP
IZI=IWK(JWIGP)
IYI=(IZI-1)/NXI0+1
IXI=IZI-NXI0*(IYI-1)
CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
1 ITI,XI(IXI),YI(IYI),ZI(IZI))
82 CONTINUE
86 JWNGP=JWNGP0+2*NNGP+1-JNGP
NGP1=IWK(JWNGP)
IF(NGP1.EQ.0) GO TO 89
JIG1MX=JIG1MN-1
JIG1MN=JIG1MN-NGP1
DO 87 JIGP=JIG1MN,JIG1MX
JWIGP=JWIGP0+JIGP
IZI=IWK(JWIGP)
IYI=(IZI-1)/NXI0+1
IXI=IZI-NXI0*(IYI-1)
CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
1 ITI,XI(IXI),YI(IYI),ZI(IZI))
87 CONTINUE
89 CONTINUE
RETURN
C ERROR EXIT
90 WRITE (LUN,2090) MD0,NCP0,NDP0,NXI0,NYI0
RETURN
C FORMAT STATEMENT FOR ERROR MESSAGE
2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S)./
1 7H MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6,
2 10X,5HNXI =,I6,10X,5HNYI =,I6/
3 35H ERROR DETECTED IN ROUTINE IDSFFT/)
END
SUBROUTINE IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
C THIS SUBROUTINE PERFORMS TRIANGULATION. IT DIVIDES THE X-Y
C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA
C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE
C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS
C CORRESPONDING TO THE BORDER LINE SEGMENTS.
C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE. POINT NUMBERS OF THE END POINTS
C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE,
C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C THIS SUBROUTINE CALLS THE IDXCHG FUNCTION.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD = ARRAY OF DIMENSION NDP CONTAINING THE
C X COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE
C Y COORDINATES OF THE DATA POINTS.
C THE OUTPUT PARAMETERS ARE
C NT = NUMBER OF TRIANGLES,
C IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE
C POINT NUMBERS OF THE VERTEXES OF THE (IT)TH
C TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND,
C (3*IT-1)ST, AND (3*IT)TH ELEMENTS,
C IT=1,2,...,NT,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = INTEGER ARRAY OF DIMENSION 6*NDP, WHERE THE
C POINT NUMBERS OF THE END POINTS OF THE (IL)TH
C BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE
C NUMBER ARE TO BE STORED AS THE (3*IL-2)ND,
C (3*IL-1)ST, AND (3*IL)TH ELEMENTS,
C IL=1,2,..., NL.
C THE OTHER PARAMETERS ARE
C IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
C INTERNALLY AS A WORK AREA,
C IWP = INTEGER ARRAY OF DIMENSION NDP USED
C INTERNALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
C WORK AREA.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),IPT(585),IPL(600),
1 IWL(1800),IWP(100),WK(100)
DIMENSION ITF(2)
DATA RATIO/1.0E-6/, NREP/100/, LUN/6/
C STATEMENT FUNCTIONS
DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
SIDE(U1,V1,U2,V2,U3,V3)=(V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
C PRELIMINARY PROCESSING
10 NDP0=NDP
NDPM1=NDP0-1
IF(NDP0.LT.4) GO TO 90
C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
20 DSQMN=DSQF(XD(1),YD(1),XD(2),YD(2))
IPMN1=1
IPMN2=2
DO 22 IP1=1,NDPM1
X1=XD(IP1)
Y1=YD(IP1)
IP1P1=IP1+1
DO 21 IP2=IP1P1,NDP0
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
IF(DSQI.EQ.0.0) GO TO 91
IF(DSQI.GE.DSQMN) GO TO 21
DSQMN=DSQI
IPMN1=IP1
IPMN2=IP2
21 CONTINUE
22 CONTINUE
DSQ12=DSQMN
XDMP=(XD(IPMN1)+XD(IPMN2))/2.0
YDMP=(YD(IPMN1)+YD(IPMN2))/2.0
C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
C NUMBERS IN THE IWP ARRAY.
30 JP1=2
DO 31 IP1=1,NDP0
IF(IP1.EQ.IPMN1.OR.IP1.EQ.IPMN2) GO TO 31
JP1=JP1+1
IWP(JP1)=IP1
WK(JP1)=DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
31 CONTINUE
DO 33 JP1=3,NDPM1
DSQMN=WK(JP1)
JPMN=JP1
DO 32 JP2=JP1,NDP0
IF(WK(JP2).GE.DSQMN) GO TO 32
DSQMN=WK(JP2)
JPMN=JP2
32 CONTINUE
ITS=IWP(JP1)
IWP(JP1)=IWP(JPMN)
IWP(JPMN)=ITS
WK(JPMN)=WK(JP1)
33 CONTINUE
C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
35 AR=DSQ12*RATIO
X1=XD(IPMN1)
Y1=YD(IPMN1)
DX21=XD(IPMN2)-X1
DY21=YD(IPMN2)-Y1
DO 36 JP=3,NDP0
IP=IWP(JP)
IF(ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21).GT.AR)
1 GO TO 37
36 CONTINUE
GO TO 92
37 IF(JP.EQ.3) GO TO 40
JPMX=JP
JP=JPMX+1
DO 38 JPC=4,JPMX
JP=JP-1
IWP(JP)=IWP(JP-1)
38 CONTINUE
IWP(3)=IP
C FORMS THE FIRST TRIANGLE. STORES POINT NUMBERS OF THE VER-
C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
C THE IPL ARRAY.
40 IP1=IPMN1
IP2=IPMN2
IP3=IWP(3)
IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
1 .GE.0.0) GO TO 41
IP1=IPMN2
IP2=IPMN1
41 NT0=1
NTT3=3
IPT(1)=IP1
IPT(2)=IP2
IPT(3)=IP3
NL0=3
NLT3=9
IPL(1)=IP1
IPL(2)=IP2
IPL(3)=1
IPL(4)=IP2
IPL(5)=IP3
IPL(6)=1
IPL(7)=IP3
IPL(8)=IP1
IPL(9)=1
C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
50 DO 79 JP1=4,NDP0
IP1=IWP(JP1)
X1=XD(IP1)
Y1=YD(IP1)
C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
IP2=IPL(1)
JPMN=1
DXMN=XD(IP2)-X1
DYMN=YD(IP2)-Y1
DSQMN=DXMN**2+DYMN**2
ARMN=DSQMN*RATIO
JPMX=1
DXMX=DXMN
DYMX=DYMN
DSQMX=DSQMN
ARMX=ARMN
DO 52 JP2=2,NL0
IP2=IPL(3*JP2-2)
DX=XD(IP2)-X1
DY=YD(IP2)-Y1
AR=DY*DXMN-DX*DYMN
IF(AR.GT.ARMN) GO TO 51
DSQI=DX**2+DY**2
IF(AR.GE.(-ARMN).AND.DSQI.GE.DSQMN) GO TO 51
JPMN=JP2
DXMN=DX
DYMN=DY
DSQMN=DSQI
ARMN=DSQMN*RATIO
51 AR=DY*DXMX-DX*DYMX
IF(AR.LT.(-ARMX)) GO TO 52
DSQI=DX**2+DY**2
IF(AR.LE.ARMX.AND.DSQI.GE.DSQMX) GO TO 52
JPMX=JP2
DXMX=DX
DYMX=DY
DSQMX=DSQI
ARMX=DSQMX*RATIO
52 CONTINUE
IF(JPMX.LT.JPMN) JPMX=JPMX+NL0
NSH=JPMN-1
IF(NSH.LE.0) GO TO 60
C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
NSHT3=NSH*3
DO 53 JP2T3=3,NSHT3,3
JP3T3=JP2T3+NLT3
IPL(JP3T3-2)=IPL(JP2T3-2)
IPL(JP3T3-1)=IPL(JP2T3-1)
IPL(JP3T3) =IPL(JP2T3)
53 CONTINUE
DO 54 JP2T3=3,NLT3,3
JP3T3=JP2T3+NSHT3
IPL(JP2T3-2)=IPL(JP3T3-2)
IPL(JP2T3-1)=IPL(JP3T3-1)
IPL(JP2T3) =IPL(JP3T3)
54 CONTINUE
JPMX=JPMX-NSH
C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
60 JWL=0
DO 64 JP2=JPMX,NL0
JP2T3=JP2*3
IPL1=IPL(JP2T3-2)
IPL2=IPL(JP2T3-1)
IT =IPL(JP2T3)
C - - ADDS A TRIANGLE TO THE IPT ARRAY.
NT0=NT0+1
NTT3=NTT3+3
IPT(NTT3-2)=IPL2
IPT(NTT3-1)=IPL1
IPT(NTT3) =IP1
C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
IF(JP2.NE.JPMX) GO TO 61
IPL(JP2T3-1)=IP1
IPL(JP2T3) =NT0
61 IF(JP2.NE.NL0) GO TO 62
NLN=JPMX+1
NLNT3=NLN*3
IPL(NLNT3-2)=IP1
IPL(NLNT3-1)=IPL(1)
IPL(NLNT3) =NT0
C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
C - - LINE SEGMENTS.
62 ITT3=IT*3
IPTI=IPT(ITT3-2)
IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63
IPTI=IPT(ITT3-1)
IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63
IPTI=IPT(ITT3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
63 IF(IDXCHG(XD,YD,IP1,IPTI,IPL1,IPL2).EQ.0) GO TO 64
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
IPT(ITT3-2)=IPTI
IPT(ITT3-1)=IPL1
IPT(ITT3) =IP1
IPT(NTT3-1)=IPTI
IF(JP2.EQ.JPMX) IPL(JP2T3)=IT
IF(JP2.EQ.NL0.AND.IPL(3).EQ.IT) IPL(3)=NT0
C - - SETS FLAGS IN THE IWL ARRAY.
JWL=JWL+4
IWL(JWL-3)=IPL1
IWL(JWL-2)=IPTI
IWL(JWL-1)=IPTI
IWL(JWL) =IPL2
64 CONTINUE
NL0=NLN
NLT3=NLNT3
NLF=JWL/2
IF(NLF.EQ.0) GO TO 79
C - IMPROVES TRIANGULATION.
70 NTT3P3=NTT3+3
DO 78 IREP=1,NREP
DO 76 ILF=1,NLF
ILFT2=ILF*2
IPL1=IWL(ILFT2-1)
IPL2=IWL(ILFT2)
C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
C - - THE FLAGGED LINE SEGMENT.
NTF=0
DO 71 ITT3R=3,NTT3,3
ITT3=NTT3P3-ITT3R
IPT1=IPT(ITT3-2)
IPT2=IPT(ITT3-1)
IPT3=IPT(ITT3)
IF(IPL1.NE.IPT1.AND.IPL1.NE.IPT2.AND.
1 IPL1.NE.IPT3) GO TO 71
IF(IPL2.NE.IPT1.AND.IPL2.NE.IPT2.AND.
1 IPL2.NE.IPT3) GO TO 71
NTF=NTF+1
ITF(NTF)=ITT3/3
IF(NTF.EQ.2) GO TO 72
71 CONTINUE
IF(NTF.LT.2) GO TO 76
C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
C - - ON THE LINE SEGMENT.
72 IT1T3=ITF(1)*3
IPTI1=IPT(IT1T3-2)
IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73
IPTI1=IPT(IT1T3-1)
IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73
IPTI1=IPT(IT1T3)
73 IT2T3=ITF(2)*3
IPTI2=IPT(IT2T3-2)
IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74
IPTI2=IPT(IT2T3-1)
IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74
IPTI2=IPT(IT2T3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
74 IF(IDXCHG(XD,YD,IPTI1,IPTI2,IPL1,IPL2).EQ.0)
1 GO TO 76
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
IPT(IT1T3-2)=IPTI1
IPT(IT1T3-1)=IPTI2
IPT(IT1T3) =IPL1
IPT(IT2T3-2)=IPTI2
IPT(IT2T3-1)=IPTI1
IPT(IT2T3) =IPL2
C - - SETS NEW FLAGS.
JWL=JWL+8
IWL(JWL-7)=IPL1
IWL(JWL-6)=IPTI1
IWL(JWL-5)=IPTI1
IWL(JWL-4)=IPL2
IWL(JWL-3)=IPL2
IWL(JWL-2)=IPTI2
IWL(JWL-1)=IPTI2
IWL(JWL) =IPL1
DO 75 JLT3=3,NLT3,3
IPLJ1=IPL(JLT3-2)
IPLJ2=IPL(JLT3-1)
IF((IPLJ1.EQ.IPL1.AND.IPLJ2.EQ.IPTI2).OR.
1 (IPLJ2.EQ.IPL1.AND.IPLJ1.EQ.IPTI2))
2 IPL(JLT3)=ITF(1)
IF((IPLJ1.EQ.IPL2.AND.IPLJ2.EQ.IPTI1).OR.
1 (IPLJ2.EQ.IPL2.AND.IPLJ1.EQ.IPTI1))
2 IPL(JLT3)=ITF(2)
75 CONTINUE
76 CONTINUE
NLFC=NLF
NLF=JWL/2
IF(NLF.EQ.NLFC) GO TO 79
C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
JWL=0
JWL1MN=(NLFC+1)*2
NLFT2=NLF*2
DO 77 JWL1=JWL1MN,NLFT2,2
JWL=JWL+2
IWL(JWL-1)=IWL(JWL1-1)
IWL(JWL) =IWL(JWL1)
77 CONTINUE
NLF=JWL/2
78 CONTINUE
79 CONTINUE
C REARRANGES THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE.
80 DO 81 ITT3=3,NTT3,3
IP1=IPT(ITT3-2)
IP2=IPT(ITT3-1)
IP3=IPT(ITT3)
IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
1 .GE.0.0) GO TO 81
IPT(ITT3-2)=IP2
IPT(ITT3-1)=IP1
81 CONTINUE
NT=NT0
NL=NL0
RETURN
C ERROR EXIT
90 WRITE (LUN,2090) NDP0
GO TO 93
91 WRITE (LUN,2091) NDP0,IP1,IP2,X1,Y1
GO TO 93
92 WRITE (LUN,2092) NDP0
93 WRITE (LUN,2093)
NT=0
RETURN
C FORMAT STATEMENTS
2090 FORMAT(1X/23H *** NDP LESS THAN 4./8H NDP =,I5)
2091 FORMAT(1X/29H *** IDENTICAL DATA POINTS./
1 8H NDP =,I5,5X,5HIP1 =,I5,5X,5HIP2 =,I5,
2 5X,4HXD =,E12.4,5X,4HYD =,E12.4)
2092 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS./
1 8H NDP =,I5)
2093 FORMAT(35H ERROR DETECTED IN ROUTINE IDTANG/)
END
FUNCTION IDXCHG(X,Y,I1,I2,I3,I4)
C THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO
C TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION
C BY C. L. LAWSON.
C THE INPUT PARAMETERS ARE
C X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA
C POINTS,
C I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1, P2,
C P3, AND P4 THAT FORM A QUADRILATERAL WITH P3
C AND P4 CONNECTED DIAGONALLY.
C THIS FUNCTION RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EX-
C CHANGE IS NECESSARY, AND 0 (ZERO) OTHERWISE.
C DECLARATION STATEMENTS
DIMENSION X(100),Y(100)
EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ),
1 (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ)
C PRELIMINARY PROCESSING
10 X1=X(I1)
Y1=Y(I1)
X2=X(I2)
Y2=Y(I2)
X3=X(I3)
Y3=Y(I3)
X4=X(I4)
Y4=Y(I4)
C CALCULATION
20 IDX=0
U3=(Y2-Y3)*(X1-X3)-(X2-X3)*(Y1-Y3)
U4=(Y1-Y4)*(X2-X4)-(X1-X4)*(Y2-Y4)
IF(U3*U4.LE.0.0) GO TO 30
U1=(Y3-Y1)*(X4-X1)-(X3-X1)*(Y4-Y1)
U2=(Y4-Y2)*(X3-X2)-(X4-X2)*(Y3-Y2)
A1SQ=(X1-X3)**2+(Y1-Y3)**2
B1SQ=(X4-X1)**2+(Y4-Y1)**2
C1SQ=(X3-X4)**2+(Y3-Y4)**2
A2SQ=(X2-X4)**2+(Y2-Y4)**2
B2SQ=(X3-X2)**2+(Y3-Y2)**2
C3SQ=(X2-X1)**2+(Y2-Y1)**2
S1SQ=U1*U1/(C1SQ*AMAX1(A1SQ,B1SQ))
S2SQ=U2*U2/(C2SQ*AMAX1(A2SQ,B2SQ))
S3SQ=U3*U3/(C3SQ*AMAX1(A3SQ,B3SQ))
S4SQ=U4*U4/(C4SQ*AMAX1(A4SQ,B4SQ))
IF(AMIN1(S1SQ,S2SQ).LT.AMIN1(S3SQ,S4SQ)) IDX=1
30 IDXCHG=IDX
RETURN
END